Bi-Laplacian Growth Patterns in Disordered Media 
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Experiments in quasi 2-dimensional geometry (Hele Shaw cells) in which a fluid is injected into a 
visco-elastic medium (foam, clay or associating- polymers) show patterns akin to fracture in brittle 
materials, very different from standard Laplacian growth patterns of viscous fingering. An analytic 
theory is lacking since a pre-requisite to describing the fracture of elastic material is the solution of 
the bi-Laplace rather than the Laplace equation. In this Letter we close this gap, offering a theory 
of bi-Laplacian growth patterns based on the method of iterated conformal maps. 



Pattern formation for two-phase flow instabilities has 
been intensely studied, both experimentally and theoreti- 
cally, for the displacement of a viscous fluid from between 
parallel plates or from a porous medium pj. In these 
cases the velocity field v(r) is well described by Darcy's 
law v(r) oc VP(r), where P(r) is the pressure. For in- 
compressible fluids V • v = 0, leading to the Laplace 
equation for the pressure, V 2 P(r) = 0, with appropri- 
ate boundary condition on the boundary of the growing 
pattern and at "infinity" . The theory for such "Lapla- 
cian Growth" patterns in 2-dimensions (i.e. r = (x,y)) 
naturally focuses on analytic functions (or their confor- 
mal inverse), simply because the Cauchy-Riemann con- 
ditions imply that the general solution of the Laplace 
equation is given by the real part of an analytic function, 
P = 'St{F(z)}, where F(z) is the unique analytic func- 
tion that satisfies the boundary conditions, and z = x+iy 

Sporadically, over the last decade, there appeared ex- 
perimental studies in which a low viscosity fluid displaces 
not a more viscous fluid, but rather a medium which is 
visco-elastic, like foam ||, clay ||, or a solution of as- 
sociating polymers ||. Elastic media are expected to 
be invaded by fracture, rather than a displacement, and 
indeed the growth patterns reported in the experiments 
had features akin to fracture patterns in brittle materi- 
als, see Fig. § Detailed comparisons with theory were 
lacking however, since an appropriate analytic theory did 
not exist. As is well known, (and see below for details), in 
fracture the relevant equation to solve is the bi-Laplace 
equation V 2 V 2 x = with appropriate boundary condi- 
tions 0. The general solution is no longer the real part 
of an analytic function, but rather 




FIG. 1: Typical pattern when water is injected into a radial 
Hele-Shaw cell filled with a solution of associated polymers 



To set up the model imagine a 2-dimensional elastic 
medium with a hole of an arbitrary shape, whose bound- 
ary z(s) is parametrized by the arc-length variable s. Into 
this hole one pushes quasistatically a fluid of pressure P. 
In equilibrium with this pressure the elastic medium will 
suffer a displacement field u(r). The strain tensor ejk 
which results is 
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(djUk + d k Uj) 



(2) 



In linear elasticity theory jjj the stress tensor is related 
to the strain tensor by 



x(z,z)=$l ztj)(z) + i/j(z) 



(1) 
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where 4>(z) and ip( z ) are a P an " °f analytic functions. 
Thus conformal techniques are not trivially applicable, 
and until recently there was no appropriate theoretical 
method to solve such equations with boundary condi- 
tions on an arbitrary ramified boundary. Even numeri- 
cal simulations were limited to lattice discretizations [^| , 
even though lattice anisotropy is a relevant perturbation 
changing the universality class of the growing patterns. 
Recent progress in the context of quasi-static fracture 
allows us to offer below an appropriate model for 
bi-Laplacian Growth patterns. 
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where E and a are material parameters. Equilibrium 
inside the elastic medium requires that 



^ dk (Jjk — for all 
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(4) 



The general solution of these equations in 2-dimensional 
is given by 

<JXX = , <Tyy = 8 2 X X , CT X y = ~d X yX , (5) 
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where the so called Airy potential x fulfills the bihar- 
monic equation 



V 2 V 2 x = 



(6) 



The solution is represented, as said, by Eq. (|l]). In order 
to develop the growth pattern we need to compute the 
tangent component of the stress tensor at the bound- 
ary of the pattern, since cracking proceeds only if this 
component exceeds a threshold ct c . To define this com- 
ponent and to state the conditions on the boundary of 
the growth pattern we use the local tangent and normal 
directions. With a being the angle between the tangent 
and the x-axis at z(s), we define derivatives with respect 
to the tangent and normal directions according to 



d t — cos(a) d x + sin(a) d y 
d n = cos(a) d y — sin(a) d x 



(7) 



The pressure P now must be balanced by the normal 
component of the stress: 



duX — a nn = — P = const. on the crack 



(8) 



Since in equilibrium no fluid slips along the boundary, 



dmX = °~tn = Ont = on the crack 



(9) 



The normal component d n nX — a tt is not determined 
from the boundary conditions, but is computedby solving 
for xi z i z )- Using the fact that 4<9 2 x/ 'dzdz = o- xx + a yy = 
&tt + o~ nn we can immediately read from Eq. (Q) , 



<j t t(z) = P + 4M[(j)'(z)] , at the boundary 



(10) 



Once we have computed the tangent component of the 
stress, we may advance the crack if Act = a u (z) — a c > 0, 
at a speed which is (on the average) proportional to Act 

Thus to compute the tangent stress and advance the 
crack we only need to determine the function <f>(z). The 
boundary conditions (||) and (ph are expressed in terms 
of cf)(z) and ijj(z) by using (Q) in (||) and (0), to derive 
d t dzX — — -P[cos(a) + isin(a)] = —Pd t z(s), where we 
identify dtz(s) as the unit vector tangent to the bound- 
ary. Rewriting this condition as dt [dgX + P z \ = we 
obtain the boundary condition on the interface p3 



(j>{z{s)) + z(s) <f>'(z(s)) + ip(z(s)) = -P z(s) + K (11) 

where ip( z ) = ^'{z) and if is a constant that can be 
chosen zero with impunity. 

The boundary conditions at infinity are obvious, since 
all stress components have to vanish as z — > oo: 

d z2 x{z,z) -> 0, d zz x{z,z) -> asz^oo. (12) 

In light of these conditions 4>(z) must have the form 
4>(z) = i(3\z + Y^j=o u -J z ~ :l w ith Pi real- The solution 
of the stress field is invariant under the transformation 



cf>^(f> + iAz + B with A real and B a complex con- 
stant. We can use this freedom to get rid of f5\ and uq, 
and write d> in the form 



3 = 1 



(13) 



Similarly from (p2|) it follows that ip has the form 



(14) 



To proceed invoke a conformal map z = $(™)(u;) that 
maps the exterior of the unit circle in the mathematical 
plane u to the exterior of the crack in the physical plane 
z, after n growth steps. The conformal map is univalent 
by construction, with a Laurent expansion 



. F W +F W/u + F { J>/u 2 + - 



(15) 



and $( ' (uj) = lu. The arclength position s in the physical 
domain is mapped by the inverse of onto a position 
on the unit circle e = exp(i0). We will be able to compute 
the stress tensor on the boundary of the crack in the 
physical domain by performing the calculation on the 
unit circle. In other words we will compute o~tt(0) on the 
unit circle in the mathematical plane. 

We perform the calculation iteratively, taking the 
stress as known for the crack after n — 1 fracture events. 
In order to implement the nth cracking event with av- 
erage velocity proportional to Act, we should choose po- 
tential positions on the interface more often when Aa(9) 
is larger. We construct a probability density P[6) on the 
unit circle e l9 which satisfies 



P(6) 



\<S>' ( - n -^(e l9 )\Aa(9)e(Aa{6)) 
|$'(n-i)( e «)|A<T(0)9(A(7(fiO)dfii 



(16) 



where 0(A<r(0)) is the Heaviside function, and 
|$ (e l9 )\ is simply the Jacobian of the transforma- 
tion from mathematical to physical plane. The next 
growth position, 9 n in the mathematical plane, is cho- 
sen randomly with respect to the probability P(9)d9. At 
the chosen position on the crack, i.e. z — <i>( ,l-1 )(e l61 "), 
we want to advance the crack with a given step of fixed 
length VAo- We achieve growth with an auxiliary con- 
formal map (j)\ n> g n (Lo) that maps the unit circle to a unit 
circle with a semi-circular bump of area A„ centered at 
e l9n |l4|, [l5|. To ensure a fixed size step in the physical 
domain we choose 



An — 



A 



'(e*»)| s 



(17) 



Finally the updated conformal map is obtained as 

= $("- 1 )(0 A „, e „H) . (18) 
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FIG. 2: Typical pattern resulting from a growth of discrete 
fracture events occuring with the probability P w B(Aa)Acr. 



The recursive dynamics can be represented as itera- 
tions of the map 4>x n ,e„{w), 



= o <fi\ 2fi2 o . . . o <f>\ ni 8 n (u) 



(19) 



Every given fracture pattern is determined completely by 
the random itinerary 

We can now represent the boundary conditions ( ^0| ) in 
terms of z(s) = <I> (,l )(e): 



0($W(e)) + *(e) ffl^M+^($(»)(e)) = -P$W(e) . 
$(«)'(e) 

(20) 

to solve this equation we introduce a power expansion for 
the ratio 



$(»)( e ) 
$(")'(e) 



E 

j=— oo 



(21) 



Note that this expansion contains both positive and neg- 
ative powers of e, whereas Eqs. (|l3| ) and ( [Ti]) contained 
only negative powers. Nevertheless upon substituting all 
the power expansions into Eq. (^0|) one finds that the 
determination of the function 4>{z) only requires the neg- 
ative powers of e (l(J , with the coefficients satisfying the 
system of equations 



U-l 



(22) 



After separating real from imaginary parts one finds an 
infinite system of linear equations. In practice we trun- 
cate at jmax = 100 and test for convergence by increasing 
the order. Note that the highest resolved Uj max requires 
computing the Fourier series (|l]) to order 2j max + 1. 

Implementing this procedure with <&(°)(o;) = to, and 
choosing Ao = Ewe generate a typical fracture pattern 
as seen in Fig. g. What is seen is the map $( 4500 )(e) 
which is topologically a circle. 

While we can guarantee that the pattern seen is indeed 
an exact bi-Laplacian pattern under the growth rules 




FIG. 3: Ff'VVAo of the pattern in Fig. § vs. n in a double 
logarithmic plot. From the slope of the least squares fit we 
estimate the dimension D — 1.4 ± 0.1. 



adopted here, we cannot guarantee that it is identical 
to any of the experimental patterns reported in ||, || [| . 
This stems from a few reasons. First, in many experi- 
ments there is a mixture of viscous and elastic phenom- 
ena, to the point that there are examples of a continuum 
of growth patterns depending on the relative importance 
of the two H ^| . In our theory we solved the bi-Laplacian 
equation after each growth step; this is relevant in the 
purely elastic limit. Secondly, and not less importantly, 
we advanced the pattern where allowed (Act > 0) at a ve- 
locity that is proportional (on the average) to Act. While 
this is accepted by a number of authors as a reasonable 
guess for the rate of growth of a fracture pattern in the 
quasi-static limit, it is by no means derived from first 
principles or universally accepted. Needless to say, in our 
procedure we can adopt any other velocity law without 
much ado, simply by changing the probability distribu- 
tion (16). We caution the reader that one does not expect 

Tol. 



the patterns to be independent of the velocity law 
Thus a more complete theory of bi-Laplacian patterns 
calls for further collaboration between experiments and 
theory to zero in on a plausible velocity law. Before doing 
so there is a limited relevance to studying carefully the 
geometric properties of the patterns obtained with this 
velocity law or another. Notwithstanding these remarks, 
we stress that the present theory offers a very convenient 
tool for assessing the fractal dimension D of the growing 
patterns. Having a univalent conformal map as in (|lq) , 
one can invoke the rigorous "1/4" theorem. This theo- 
rem states that if R n is the radius of the minimal circle 
that contains the pattern after n fracture events, then 



(n) 



> \n n . 



(23) 



Accordingly one expects that for large n the first Laurent 
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coefficient satisfies 

F[ n) « n 1 ' . (24) 

One advantage of the present approach is that the first 
Laurent coefficient is known exactly as 

F[ n) = n n k=1 ^iTT k , (25) 

which is computable to machine precision. In Fig. ^ 
we present, in double logarithmic plot, F^ /\Ao vs. n. 
Reading the slope of the least-squares fit yields a dimen- 
sion D = 1.4 ± 0.1. 

The method of iterated conformal maps offers conver- 
gent calculations of fractal and multifractal properties of 
the growth patterns. If required, one can use the formal- 
ism to obtain highly accurate values of the fractal dimen- 
sion, cf. Era. In addition, if the properties of the growth 



probabilities are of interest from the multifractal point 
of view, there are available methods to compute these in 
a convergent scheme Il7j that is not available in direct 
numerical simulations. However, such refinements would 
be justified only after future work to solidify further the 
relation between theory and experiment. 
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